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ABSTRACT 



To date, there is no core accretion simulation that can successfully account 
for the formation of Uranus or Neptune within the observed 2-3 Myr lifetimes 
of protoplanetary disks. Since solid accretion rate is directly proportional to the 
available planetesimal surface density, one way to speed up planet formation is to 
take a full accounting of all the planetesimal-forming solids present in the solar 
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nebula. By combining a viscously evolving protostellar disk with a kinetic model 
of ice formation, which includes not just water but methane, ammonia, CO and 
54 minor ices, we calculate the solid surface density of a possible giant planet- 
forming solar nebula as a function of heliocentric distance and time. Our results 
can be used to provide the starting planetesimal surface density and evolving solar 
nebula conditions for core accretion simulations, or to predict the composition of 
planetesimals as a function of radius. 

We find three effects that favor giant planet formation by the core accretion 
mechanism: (1) a decretion flow that brings mass from the inner solar nebula to 
the giant planet-forming region, (2) the fact that the ammonia and water ice lines 



should coincide, according to recent lab results from Collings et al. (2004), and 
(3) the presence of a substantial amount of methane ice in the trans-Saturnian 
region. Our results show higher solid surface densities than assumed in the core 
accretion models of Pollack et al. (1996) by a factor of 3-4 throughout the trans- 
Saturnian region. We also discuss the location of ice lines and their movement 
through the solar nebula, and provide new constraints on the possible initial disk 
configurations from gravitational stability arguments. 

Keywords: solar nebula; planetary formation; ices 



1. Introduction 

Previous core accretion simulations of the formation of Uranus and Neptune require the 
solar nebula lifetime to be at least 10 Myr longer than the observed median protostellar 
disk lifetime (2-3 Myr; Haisch et al. 2001). Similarly, the simulation of Saturn's formation 
by Alibert et al. (2001) allows the planet to form in 2.5 Myr, but only if it migrates 
from 11.9 AU to 9.5 AU during formation. This scenario is incompatible with the outward 
migration of Saturn predicted by the Nice model of planetary dynamics (Tsiganis et al. 2005) 
This research was motivated by the need to create a model of the solar nebula with enough 
solid mass to form Jupiter, Saturn, Uranus and Neptune via core accretion at heliocentric 
distances of 5 < R < 20 AU. 

For core accretion simulations, planetesimal mass available for building giant planet 
cores is usually calculated by invoking the scaled minimum-mass solar nebula with surface 
density £ oc i?~ 3 / 2 , obtained from spreading out the current mass of the planets (Wei- 
denschilling 1977), and assuming a gas/solid ratio of ~ 70. However, a single planet core 
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embedded in a disk of planetesimals grows at the rate of 



core 



7ra core^solid^-^c 



(1) 



dt 



where F g is a factor taking into account gravitational focusing by the core (Safronov 1969), 
Ssoiid is the solid surface density of planetesimals, fl is the orbital angular frequency of the 
core, and a core is the core's effective capture radius. For Keplerian orbits, we have 



The i?™ 3 / 2 dependence of f2 means that even planets forming in a disk with uniform surface 
density should, in the end, reach masses that suggest £ oc R~ 3 ^ 2 if growth rate is not taken 
into account. A near-uniform surface density distribution, in which mass is not concentrated 
in the inner disk, gives higher surface densities in the outer disk and may speed up the 
formation of Uranus and Neptune. 

Furthermore, the canonical gas/solid ratio of 70, which is based on the composition of 
Comet Halley (Jessberger et al. 1989), includes only one ice in the solid inventory — H 2 0. 
However, there were certainly other ices present in the solar nebula: the Deep Impact ejecta 
from Tempel 1 contained ~ 1000 tons of excavated CO, CO2 and CH3OH (A'Hearn 2008). 
Taking a full accounting of the solids in the solar nebula may reveal other volatile mass 
sources besides H 2 that condense in the outer solar nebula and improve the formation 
prospects of Uranus and Neptune. 

In this work, we calculate the time-evolving solid surface density available for giant 
planet formation. We make two updates to previous solar nebula simulations: (1) we relax 
the assumption of a steady-state disk and allow for viscous redistribution of mass, which 
tends to drive the surface density profile toward uniformity, and (2) we add a chemical 
reaction network to trace the formation and freezeout of volatiles. 

This paper is organized as follows. In £[2j we outline our method for calculating solid 
surface density in the solar nebula: which codes we use, how they are linked, and how they 
differ from previous solar nebula models. In £[3j we describe our protostellar disk model and 
use it to calculate the total surface density (gas+solid) and temperature evolution of the 
disk. §[1] details the physics behind our chemical reaction network, the abundances of the 
major ice species, and the location of the ice lines. We calculate the solid surface density 
available for planet formation in £j5] and present our conclusions in §[6} 




(2) 
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2. Simulation Overview 

This simulation consists of two independent, non-interacting codes used in sequence. 
The first code is a non-steady-state, 1+1-dimensional a-disk model that simulates the dy- 
namical evolution of the solar nebula. Code 1 gives the surface density distribution t), 
along with the vertical and radial temperature and volume density distributions T(R,z,t) 
and p(R, z, t). 

The results from Code 1 fill an important niche in the body of work on protostellar disk 
evolution. Currently available are 2-d, static disk models that include stellar insolation and 
radiogenic heating (e.g. D'Alessio et al. 2006); quasi-steady-state 1+1-d models where M(t) 
follows an imposed prescription (e.g. Hersant et al. 2001); and non-steady-state models with 
only one zone in the vertical direction (Ruden and Pollack 1991). The evolving protostellar 
disks of Boss (2007) include 3-d hydrodynamics and radiative diffusion, but cover a timespan 
of < 4000 years. 

In this work, we couple full vertical structure models with the time- dependent, non- 
steady-state surface density evolution of the disk and present the first long-term solar neb- 
ula simulation with viscous angular momentum transport determined self-consistently with 
vertical structure. For a description of how the resulting dynamical evolution of our disk 



differs from previous simulations, see £3.7.2 



Code 2 is a non-equilibrium chemical model, including gas reactions, grain-surface re- 
actions, freezeout and desorption. This code traces the abundance of 211 species, including 
58 ices, as the disk cools and the ice lines move inward. Chemical reaction, freezeout and 
desorption rates are governed by the temperature and density of the solar nebula gas and 
grains. We use the midplane temperatures and densities calculated by Code 1 to generate 
reaction rate coefficients for Code 2. 

Previous chemical models of the solar nebula include those of Aikawa et al. (1996), 



Willacy and Langer (2000), and Ilgner and Nelson (2006). Aikawa et al. (1996) extended 
the minimum-mass solar nebula (MMSN; Hayashi 1981) to 800 AU to predict the observed 
depletion of gaseous CO in the outer regions of T-Tauri disks. Their calculation included 
gas-phase reactions of CO and other C, H, O, N and S-based molecules, freezeout of CO and 
its reaction products, formation of H 2 , and electron-ion recombination on grain surfaces. 
The density and temperature of the underlying disk were assumed to stay constant with 
time, while the time-evolving midplane CO abundance was calculated with an independent 
sequence of radial one-zone models. 

Willacy and Langer (2000) calculated the vertical distribution and column density of 
common interstellar molecules as a function of time in the flared protoplanetary disk of 
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Chiang and Goldreich (1997). They used a similar species set to Aikawa et al. (1996) 
but added grain-surface reactions and photodesorption to the chemistry. The simulation 
consisted of a 2-d grid of noninteracting one-zone models in the R — z plane. The underlying 
temperature and density distribution of the protostellar disk remained constant with time. 

Finally, Ilgner and Nelson (2006) embedded a five-species chemical reaction network in 
an evolving a-disk model, fully coupling the chemical and dynamical evolution of the disk. 
They used the chemical reaction network to trace the ionization fraction, which is affected by 
charge-transfer reactions, electron-ion recombination, and viscous diffusion. The simulation 
was computationally tractable due to the small reaction network (5 species, 4 reactions), 
and the limited extent of the protostellar disk (0.1-10 AU). 

Code 2 of this work is the first to combine an extensive chemical reaction network, 
including ice formation and grain-surface reactions, with a time-evolving protostellar disk 
model (as used by Ilgner and Nelson 2006). Our chemical model consists of a radial sequence 
of noninteracting one-zone models, as in Aikawa et al. (1996). The disk model in Code 1 
therefore functions as a stand-alone input to the chemical model of Code 2. Our models 
differ from those of Aikawa et al. (1996) and Willacy and Langer (2000) in that we allow 
the disk temperature and density to change with time, as prescribed by the results of Code 
1. Although Code 1 simulates the full vertical structure of the protostellar disk, we apply 
Code 2 only to the midplane, which is where solids settle and build planetesimals. 

Code 2 calculates the abundance of each ice and gas species as a function of heliocentric 
radius and time. Combining these results with information from the literature about the 
composition of refractory grains, we calculate the gas/solid mass ratio in the solar nebula. 
This ratio tells us what fraction of the total surface density calculated by Code 1 is solid 
and can build giant planet cores. The end result of both codes is the solid surface density 
distribution E so nd(-R, t). 

Since most molecules in the solar nebula formed in a molecular cloud, we first run Code 
2 for molecular cloud conditions (T = 10 K, hydrogen number density 2 x 10 4 cm -3 , grain 
radius 0.258 fim). The molecular cloud model begins with atomic gas of solar composition 
(except for hydrogen, 99% of which is in H2) and bare refractory grains at the standard 
interstellar abundance n g = 10~ 12 (nn + 2nn2), where n g is the grain abundance and nu and 
nH2 are the abundances of H and H 2 , respectively. We run the molecular cloud model for 
10 6 yr and use the resulting ice-gas mixture as the starting composition of the solar nebula. 
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3. Disk Model 

3.1. Overall Model Formulation 

The goals of this simulation are to calculate the time-evolving temperature and density 
in the solar nebula, which determine the reaction rates for our kinetic model of ice formation, 
and to follow the viscous evolution of the surface density profile. Although from a planet 
formation perspective we are primarily concerned with the ice inventory at the solar nebula 
midplane, an accurate dynamical description of the midplane requires understanding the 
vertical energy balance in the disk. To simplify the problem of long-term protostellar disk 
evolution, we use the following assumptions: 

1. The disk is axisymmetric and symmetric about the midplane. 

2. The disk is geometrically thin: H/R « 1. 

The disk is represented in cylindrical polar coordinates R (heliocentric distance) and z 
(height above midplane), where the z-axis is the rotation axis of the nebula. Assumption 
1 reduces the physical 3-dimensional disk to a 2-dimensional quadrant with zero flux at 
the midplane, allowing us to suppress the azimuthal coordinate 9. Assumption 2 allows 
the vertical and radial dimensions of the disk to be decoupled to form a 1+1-dimensional 
framework: at each radial gridpoint resides an independent vertical structure model. At 
each timestep, energy balance between viscous heating and radiative cooling is solved within 
in a single annulus. 



3.2. Radial Diffusion 



We calculate the radial motion of mass using the surface density diffusion equation of 
Lynden-Bell and Pringle (1974): 



~dt 



rM 



Rl/2 M 



(3) 



Eq. [3]is valid when the protostar's accretion rate is small, so the sun is almost fully assembled 
at the beginning of the simulation. This assumption is consistent with the T-Tauri phase of 
disk evolution. We take the viscosity v governing radial diffusion to be the midplane value. 



To calculate f m id p i a 
transport in the disk: 



we make the following assumptions about angular momentum 
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1. Viscous stresses follow the standard a-viscosity prescription (Shakura and Syunyaev 
1973). We assume the magnetorotational instability (Balbus and Hawley 1991; here- 
after MRI) provides the requisite turbulent viscosity. 

2. The grains and gas are well mixed. This assumption allows us to track the disk's solid 
surface density evolution by modeling only a single fluid. 



The turbulent viscosity is 



v = -ac s H, 



(4) 



where c s is the isothermal sound speed and H is the modified pressure scale height, softened 
into a nonsingular form (Milsom et al. 1994): 



H 



c s /Q 



y/1 + (2*W/<*) 

The angular velocity, neglecting disk self-gravity, is 

1/2 



(5) 



GM* 



and the isothermal sound speed is 



(R 2 + z 



Rg 



2^,3/2 



T. 



(6) 



(7) 



In parts of the disk where H > z sur{ , the disk surface height, we substitute z suri for H in Eq. 

We take a mean molecular weight of fi = 2.33 g mol" 1 (Ruden and Pollack 1991). Based 
on recent global-disk MHD simulations (Lyra et al. 2008), we choose a = 2 x 10~ 3 . We 
further discuss the interplay between a, disk mass and surface density profile, along with 



our constraints on these parameters, in £3.6. All free parameters in the disk model are listed 
in Table [TJ 

Although we assume turbulence is caused by MRI, we do not include the radial variation 
of the viscosity coefficient a caused by the dead zone (e.g. Kretke and Lin 2007, Reyes-Ruiz 
et al. 2003). However, we do account for the cessation of MRI turbulence in the tenuous 



outer regions of the disk, where ions and neutrals may be only weakly coupled. Hawley 



and Stone (1998) find that neutrals participate fully in MRI turbulence only when the ion- 



neutral collision timescale T co u is much less than the orbital timescale; T co u < 0.01/fi. Ions 
and neutrals decouple and MRI turbulence ceases when T coll > 100/ ft. In the transition 
region between coupled and decoupled neutral-ion interaction, where 0.01 < T co nQ < 100, 
we decrease log a linearly from its fiducial value. 
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To calculate T co ii, we find the mean free time of a neutral particle traveling through a 
gas of moving ions. The neutrals and ions collide at an average angle of 90° and thus a mean 
impact velocity of 



«coU = \l Vl + v\ = V n I 1 + 



1/2 



where fi n and ft are the molecular weights of neutrals and ions, respectively. The mean 
collision time is then 

Tco11 = ~i , 7 777 / 772' ( 9 ) 

7T (a n + ft) UiV n (1 + fin/ ft) 

where a n and a, are the radii of neutrals and ions, v n = V3c s is the mean velocity of a 
neutral and ft is the number density of ions: 

n t = I* (10) 

We assume an ionization fraction of fi = 10~ 9 and an ion mass of m; = 30 AMU. fi and rrii 
are listed along with other free parameters in the simulation in Table [TJ 



3.3. Vertical Structure 

Since the disk midplane viscosity is determined by the balance between viscous energy 
generation and radiative energy losses at the photosphere, we calculate full vertical structure 
models at each radial grid point. We use the following assumptions about the mass and 
energy balance in the vertical direction: 

1. The disk stays in hydrostatic equilibrium. Even though the mass distribution evolves, 
the vertical sound crossing timescale is much less than the radial diffusion timescale. 

2. We neglect stellar irradiation as an energy source, as it has little effect on the midplane 



temperature. In section 3.6 we show that the disk is self-shadowed (H/R decreases with 



R), making this assumption self-consistent with our simulation results. 

A system of three coupled differential equations specifies the vertical structure of the 
disk. First, the disk is in hydrostatic equilibrium: 

f z - -fit*. (") 

Second, the viscous energy generation rate per unit volume (Pringle 1981) is 

BF 9 
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Third, where vertical energy transport is regulated by radiative diffusion and energy trans- 
port in the radial direction is negligible (thin disk assumption), the energy equation is 

dT OPT , , 

& = a7 v F < 13 > 

To calculate the thermodynamic gradient V = d In T/d In P, we use the Schwarzschild 
criterion for stability against convection (Kippenhahn and Weigert 1994): 

y _ j V r ad, V ra d < V a d 
V ra d > V a d 



In Eq. 14, V a d is the adiabatic gradient, which is 2/7 for a diatomic gas; V con v is the 
convective gradient; and V ra d is the radiative gradient: 

3 kPF 



Vrad "4ac^T 4 ' (15) 

where a is the radiation density constant and k is the local Rosseland mean opacity (for a 
description of the opacity tables used in this calculation, see £3.5). 

We use mixing-length theory to calculate V con v We begin the convection treatment by 
defining two dimensionless quantities: 

W EE V rad - Vad (16) 

and 

3acT 3 r— 

U = M^- (17) 

In Equation 17, cp is the specific heat at constant pressure and l m is the mixing length. 

As in Milsom et al. (1994), we define i m as the minimum of the convection zone top or 

the pressure scale height H (Equation [5]). e is an analytical softening of the ratio H/Q 2 z, 
eliminating the singularity at z = 0: 



n 4 (l + z 2 /H 2 )' 



Defining a new variable £ as 



£ EE V conv - V ad + U 2 = VVconv ~ V e + U (19) 

(where V e is the thermodynamic gradient of an individual convective element), we can then 
solve the cubic equation 

(Z-U) 3 + 8 ^-(e-U 2 -W)=0 (20) 
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to find V conv . For the derivations of Equations [T9] and [20] see Chapter 7 of |Kippenhahn and 



Weigert (1994). 



Finally, we close the system of equations with the ideal gas equation of state, 



Rg 



pT. 



(21) 



3.4. Computational Methods and Boundary Conditions 

We begin the simulation with a specified surface density profile for the disk, S oc r~ 3 / 2 
(Hayashi 1981). The total disk mass used to normalize S(r) and the viscosity coefficient 
a are determined by constraints of gravitational stability and the need for enough mass to 
form giant planets. For full details of how the t — disk is constructed, see §3.6 We also 



assume the disk is embedded in a medium with an ambient temperature T am b (for example, 
a molecular cloud core). This places a firm lower limit on the nebula temperature and forces 
inactive regions of the disk to assume an isothermal vertical structure. 



As in Bell et al. (1997), we place the disk surface at optical depth r = 0.03. We begin 
the vertical structure calculation with guesses for the free parameters T sur f and p sur f, the 
density and temperature at the r = 0.03 surface. 

The net flux leaving the disk surface is determined by the accretion temperature T acc , 
generated solely by viscous stresses: 

F = vTL. (22) 



We require a value of T acc at the disk surface to to start integrating Equations [12] and [13} 
The free parameter T sur f is the flux sum of the accretion temperature T acc and the ambient 
temperature T am b, invoking the Eddington approximation: 

TLf = lib + \ [r + f (r)] T a 4 cc . (23) 



In equation 23, /(r) is the Hopf function, with a value of 0.601242 for r = 0.03 (Bell et al. 



1997; Mihalas 1978). 



After calculating P from Eq. 21, the variables of integration T acc , P and F are all 
specified. We find the height of the r = 0.03 surface by solving the following equation for z: 

1 _ r 

W~z ~ k( P ,T)P' ^ 

Integration proceeds inward from the disk surface to the midplane using the Runge-Kutta 
method with adaptive stepsize control (Press et al. 1992). We repeat the integration, 
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adjusting T surf and p sur f with the Newton- Raphson algorithm (Press et al. 1992), until a 
solution is found that satisfies 

-^midplane = (25) 



(following Assumption 1, £3.1) and 



2 pdz = E. (26) 

Jz=0 

After the vertical structure is solved for each radius, we use the resulting z/ mi d p i a nc to 
update S(r). We implement implicit finite differencing to solve Eq. |3j adjusting the timestep 
At so that surface density varies by a maximum of 1.5% during a single timestep (e.g. Ruden 
and Pollack 1991). 

The radial grid is equispaced in r 1 / 2 and contains 600 cells between the inner boundary 
at 0.3 AU and the initial outer boundary at 30 AU. We implement a zero viscous-stress 
inner boundary condition, so matter falls directly onto the star from the innermost gridpoint. 



Although the actual accretion shock is at a much smaller radius than 0.3 AU, Ruden and Lin 



( 1986 ) find that the location of the inner boundary R in is unimportant as long as R OVLt >> Rin- 
We allow the disk to expand freely from the original outer boundary, so new radial grid cells 
are added to the disk as the simulation progresses. 



3.5. Opacities 



Solving Eq. 15 requires knowing local mean opacity. Semenov et al. (2003) calculated 
Rosseland mean opacities due to icy grains for temperatures down to 5K. The authors 
modeled several different grain types, including spheres, aggregates, and porous grains. We 
use the 5-layered sphere topology, where each grain consists of a silicate/iron nucleus covered 
with successive layers of volatiles as the temperature decreases. Following Semenov et al., 
we assume the grains condensed from gas of solar composition, with the abundances given 
by Helling et al. (2000). These abundances are primarily derived from the catalog of Anders 
and Grevesse (1989). We list the abundances of the most common elements in Table [2j along 
with the types of solids they can form: rock, metal, ice, or refractory CHON (e.g. graphite 
or kerogen; see Jessberger et al. 1988). 

At temperatures above 1200K, iron and silicates sublimate and the main opacity source 
switches from dust to molecular gas. Semenov et al. (2003) caution that the gas opacities in 
their model are only approximate in the temperature range where dust and gas opacity are 
comparable. For temperatures above 1000K, we use the opacities of Ferguson et al. (2005), 
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updated for our chosen set of solar abundances (J. Ferguson, private communication). It 
should be noted that the number of atomic and molecular lines included in the calculation 
of an opacity table has a greater impact on the resulting opacities than minor changes in 
assumed solar composition. 

Regions where opacity is a strong function of temperature can cause convergence prob- 
lems for the Newton-Raphson algorithm, which will oscillate between two solutions with 
similar temperatures but quite different opacities. To mitigate this problem, we smoothed 
the raw opacity tables with an averaging filter. Figure [T] shows opacity as a function of 
temperature and density in our model. 



3.6. Initial Conditions 

Constraints on the early stages of the solar nebula come from observations of Class 1 
protostars and T-Tauri stars and the presence and orbital configuration of the planets. We 
require our t — disk to meet two criteria: 

1. The disk is gravitationally stable to axisymmetric perturbations at all radii. 

2. The disk is massive enough to produce Jupiter and Saturn by core accretion within 
the observed mean protostellar disk lifetime of 2-3 Myr (Haisch et al. 2001). 

To test whether a particular disk meets Criterion 1, we calculate the Toomre Q param- 
eter at every radial gridpoint. Gravitational stability requires that 

Gravitationally unstable disks will fragment and form stellar or substellar clumps on a 10 3 
year timescale. The lack of a binary companion for the sun and the near-circular, non-chaotic 
orbits displayed by all the planets makes gravitational fragmentation anytime during solar 
nebula evolution unlikely. 

Our Q > 1 criterion determines the combinations of viscosity parameter a, mass dis- 
tribution and total disk mass that create viable accretion disks. High values of a act to 
stabilize the disk by increasing energy generation and therefore sound speed. Large disk 
masses destabilize the disk by increasing E. Surface density profiles that decline steeply 
with radius damp gravitational instability by loading the mass preferentially toward the 
inner disk, where dissipative shear stresses are strongest. 
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Hersant et al. (2001; hereafter HGH) used an evolving a-disk model and measured 
solar-system D/H ratios to constrain the structure of the solar nebula. They chose the 
surface density profile £ oc R . Even with an extremely high turbulent energy generation 
rate, a > 0.01, disks with this surface density distribution are gravitationally stable only if 
Mdisk < O.O8M . However, Criterion 3 requires that there be sufficient solid mass to form 
Jupiter and Saturn within ~ 3 Myr. Assuming a gas/solid mass ratio of 70 (Pollack et al. 
1996), a disk with E oc R^ 1 out to 30 AU and M disk = O.O8M has a solid surface density 



of only 6 g cm -2 in Saturn's feeding zone at 9.5 AU. Hubickyj et al. (2005) found that a 
planet forming from 6 g cm -2 of solids at 5 AU takes more than 13 Myr to initiate rapid 
gas accretion, necessary for forming a massive atmosphere. Since planet formation proceeds 
more slowly at larger heliocentric distances (Eq. [I]), this disk is unlikely to be able to form 
Saturn. We thus require a more massive disk, with a steeper surface density profile to keep 
the outer regions gravitationally stable. 

We choose a disk mass Mdisk = O.12M and a surface density E oc R~ 3 ^ 2 (Hayashi 
1981). We take a = 2 x 10~ 3 , the approximate azimuthally and vertically averaged value 
at Jupiter's heliocentric distance in the global MHD turbulence simulations of Lyra et al. 
(2008). Our solar nebula has inner and outer radii of R{ Q = 0.3 AU and i? ou t = 30 AU, 
covering the range of heliocentric distances where planets exist today. It is worth noting 
that studies of the relative motions of dust and gas predict a much larger outer radius: the 
gas drag on cm-sized bodies can cause them to be quickly lost to the Sun unless they have 
a very large distance to migrate (Stepinski and Valageas [1997]; Ciesla and Cuzzi [2006]). 
Mitigating the deleterious effect of gas drag on large grains may require that the nebula has 
an additional low-density component outside R = 30 AU, which we do not model in this 
work. 

We assume our t = disk coincides with the beginning of the sun's T-Tauri phase. 
Most of the sun's mass is already in place: M* = O.95M at the start of our simulation. We 
assume an ambient temperature of 20 K, approximating a disk embedded in the remnants 
of a molecular cloud. All free parameters in the disk model are listed in Table [T] 

Although our model neglects solar irradiation, we find that the giant planet-forming 
region of our t=0 disk is in fact shadowed. The maximum aspect ratio H/R of the r = 2/3 
photosphere occurs at R = 2 AU. Beyond this radius, no part of the disk receives direct 
sunlight. The shadow persists until t ~ 1 Myr, when regions with R > 20 AU attain an 
aspect ratio greater than the local maximum at 2 AU. In this way, our disk evolution is 
broadly consistent with model T-Tauri disk SEDs, which require flaring to explain observed 
spectra (e.g. Chiang et al. 2001), yet the giant planet-forming region (~ 5-20 AU) still 
receives no solar heating during the icy planetesimal formation epoch. 
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3.7. Disk Model Results 

3.7.1. Surface Density 

We run the solar nebula model for 2 Myr, about the median age of observed protostellar 
disks (Haisch et al. 2001). Figure [2] shows the evolution of the disk surface density profile in 
the giant planet-forming region, 5-20 AU (Tsiganis et al. 2005). The surface density profile 
at t = is shown on the upper left. Subsequent panels show snapshots of S(-R) as the 
simulation progresses plotted in black. For reference, surface density profiles from preceding 
plot panels are retained in gray. Figure [3] shows the surface density evolution of the entire 
disk as it expands from 30 AU at t = to > 80 AU after 2 Myr. 

From examining Fig. [2j we can divide the evolution of the solar nebula's giant planet- 
forming region into two epochs. In Epoch 1, the surface density profile flattens as mass is 
redistributed from the inner disk to the outer disk. At t — 0, the net mass flow is outward 
everywhere except at R < 1 AU, which is accreting inward toward the star. A wave of 
local surface density enhancement consequently propagates outward, reaching ~ 10 AU at 
t = 3 x 10 4 yr and dissipating near ~ 17 AU at t = 3 x 10 5 yr. At the end of Epoch 1, at 
t~5x 10 5 yr, the solar nebula has assumed a surface density profile flatter than £ oc 
(equal mass in each annulus) in all but the inner 0.5 AU of the disk. The effect of this mass 
redistribution epoch, therefore, is to drive the surface density distribution toward uniformity. 

In Epoch 2, t > 5 x 10 5 yr, surface density decreases with time everywhere in the disk. 
In Fig. [3j we see the solar nebula grow more tenuous as it expands: after 2 Myr, the total 
disk mass is 0.067 M Q spread over 80 AU, as opposed to 0.12 M & within 30 AU at the 
beginning of the simulation. The protosun reaches a final mass of 1.003 M . 

Figure [4] shows surface density as a function of time at 5 AU, 10 AU, 15 AU and 20 AU. 
At each heliocentric distance, we see the temporary mass growth as the density wave passes 
through. The largest relative surface density increase AS/So occurs at 11.5 AU, where the 
initial surface density of 442 g cm -2 grows by 28% to 558 g cm -2 in 85,000 yr. This surface 
density enhancement moving outward through the giant planet-forming region strongly favors 
solid planet-core formation. 

At the beginning of the sun's T-Tauri phase, before substantial grain growth occurs, 
gas and grains are well mixed: accreting gas drags grains with it. If gas moves inward with 
a constant mass flux, surface density decreases with time everywhere in the disk and solids 
are lost into the sun. Our net outward gas flow deposits small grains in the 5-15 AU region 
within the planetesimal formation timescale of ~ 5 x 10 4 yr (Hubbard and Blackman 2006). 
These grains provide surface area for ice deposition and add solid mass to the giant planet 
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feeding zones, speeding up solid core formation. 

The existence of Epoch 1 and consequent driving of the surface density distribution 
toward uniformity is a robust result that we see even using t — surface density profiles 
with power-law indices of 0.7 and 1. Our chosen fiducial disk remains stable to axisymmetric 
perturbations (c.f. Eq. 27) throughout its evolution even as mass flows outward. 



3.7.2. Mass Flux 



In their study of the similarity solutions of viscously evolving disks, Lynden-Bell and 



Pringle (1974) found that the natural behavior of an accretion disk is to diffuse and spread 



out, creating a net inward flow in the inner disk, a turnaround radius, and an outward 



decretion flow in the outer disk. Ruden and Lin (1986) uncovered the same behavior in 



numerical simulations of the evolution of the minimum-mass solar nebula. Why haven't 
many other a-disk models predicted the outward mass flow into the trans-Saturnian region 
that could so improve the prospects for giant planet formation? Most previous evolving 



disk models such as HGH and the pseudoevolutionary sequence of Bell et al. (1997) have 



required the disk to be in steady state. In steady-state models, each annulus receives an 
equal inward-driven mass flux with the sun serving as the sink particle. There is therefore 
no chance for any part of the disk interior to -R out at t — to gain mass. Note that the 
assumption of steady-state, net inward accretion is still consistent with viscous expansion of 
the disk from its original outer boundary, as occurs in the HGH models and this work. 



The non-steady-state evolving disk model of Ruden and Pollack (1991) assumed tur- 



bulent stresses were caused solely by convection, and set the viscosity to zero in optically 
thin regions where radiative transport is efficient. Far from the star, the shallow potential 
well does not allow the disk to attain high volume densities, so reaching an optical depth 
of r > 1 requires more mass than in the inner disk. At R = 30 AU, the outer radius at 
the beginning of the Ruden and Pollack calculation, convection ceases almost immediately 
and the disk becomes quiescent. The inward accretion flow then removes enough mass from 
progressively interior annuli to make them optically thin, leaving behind a remnant surface 
density profile where £ increases with R. 

In the Ruden and Pollack simulation, mass is not redistributed from the inner disk 
to the outer disk — the slope of S(-R) changes from negative to positive because the inner 
parts of the disk remain optically thick for the longest, participate in convection-driven 
accretion the longest, and therefore lose the most mass. However, in their analysis of the 
coupling between the kinetic stress tensor governing hydrodynamic stability, the magnetic 
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stress tensor, and the shear field, Hawley et al. (1999) find that MRI is uniquely capable 



of providing turbulent angular momentum transport in disks. By relaxing the steady-state 
condition, but still permitting angular momentum transport by MRI in optically thin regions 



(though MRI ceases if ions and neutrals decouple; see £3.2), we allow outward flows to develop 
and viscous stress to redistribute mass from the inner disk to the giant planet-forming region. 

Figure [5] shows time snapshots of the mass accretion rate as a function of radius. The 
location with M = is the turnaround between inward accretion and outward "decretion." 
This turnaround radius moves outward over time, so that more and more of the disk is 
accreting toward the star. This behavior is required to open a hole in the inner disk and 
make the transition from classical T-Tauri to non-accreting, weak-lined T-Tauri star. Instead 
of accreting onto the star, the outwardly decreting material (R > 15 AU) will eventually be 
lost to photoevaporation (Alexander et al. 2006). 



3.7.3. Midplane Temperature 



Figure [6] shows shapshots of the temperature of the protostellar disk midplane at R < 
30 AU as the disk evolves. Outside of 30 AU, the entire disk (midplane and above) is 
isothermal at the ambient temperature throughout the simulation (20 K, plotted with a 
dash-dotted line). The first three panels of Fig. [6] show that during Epoch 1, areas of the 
disk that increase in mass are correspondingly heated: the temperature at 10 AU increases 
from 40 K at t = to 49 K at 3 x 10 4 yr. 

We first evaluate the suitability of our model disk for delivering crystalline silicates 



to comets, as detected by Keller et al. (2006) in comet 81P/Wild 2 dust captured by the 



Stardust spacecraft. The dashed line in Fig. [6] shows the minimum temperature for silicate 
crystallization of 800 K (Gail 1998). (Note that forsterite formation requires T > 1067 K; 
Hallenbeck et al. [2000].) The inner nebula retains temperatures above 800 K for 10 5 yr. 



Dullemond et al. (2006) calculated the formation and transport rates of crystalline silicates 



in a 1-d evolving a-disk model with infall from a collapsing cloud. Relaxing the steady- 
state assumption, they found that viscous spreading of the disk, rather than turbulent radial 
mixing, could relocate crystalline silicates produced early on in the inner, hot (T > 800 K) 
parts of the disk to distances of R > 15 AU, where comets formed. A gradual outward 
push of crystalline silicates is consistent with the nebula flow pattern observed in this work, 
and may explain the anticorrelation between disk age and crystallinity in the inner disk 



(R < 10 AU) discussed by van Boekel et al. (2004). However, the danger that crystalline 



silicates transported outward through the midplane would be swept up by protoplanets before 



ever reaching the comet-forming region led Ciesla (2007) to favor more rapid transport than 
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the decretion flow in this model. 

The dotted line on Fig. [6] shows the H 2 ice line at 160 K (Collings et al. 2004). 
Although we require a chemical kinetic model to determine the exact position and width of 
the ice line (since ice deposition is not instantaneous but proceeds at a finite rate limited by 
the gas velocity of H 2 and the available grain surface area), we can nevertheless estimate 
its location and compare with previous disk models. In our simulation, the H 2 ice line 
is always inside the formation zones of Saturn, Uranus and Neptune, and crosses Jupiter's 
feeding zone (5.2 AU) at 5 x 10 4 yr. All giant planet cores, therefore, should have formed 
from icy grains. 

After 2 Myr, the midplane outside 15 AU has become isothermal at the ambient tem- 
perature. This temperature, therefore, has a profound effect on the ice inventory of the solar 
nebula. The most volatile ices (CO, CH4 and N 2 ) have sublimation temperatures below 
45 K. Simply raising the ambient temperature above this value would ensure that none of 
these ices could freeze on grains. However, C/H enrichments of 20-40 over the solar value 
measured in Uranus and Neptune's atmospheres indicate that direct gas accretion from the 
nebula cannot explain these planets' composition (Baines et al. 1995). We therefore postu- 
late that the outer solar nebula midplane contained CO and CH4 ice, making our choice of 
ambient temperature appropriate. Further evidence for a cold solar nebula comes from the 
observed noble gas enrichment in Jupiter's atmosphere — trapping argon in icy planetesimals 
requires temperatures of 30 K (Owen et al. 1999). 

In the next section, we describe our chemical model of the solar nebula midplane. We 
will use these results, along with the surface density curves shown in Fig. [2j to predict the 
solar nebula's solid mass distribution in £j5j 

4. Chemical Model 
4.1. Reaction Set 

The foundation of our chemical model is the UMIST database RATE95 (Millar et al. 
1997). We follow the chemistry of 211 species, 153 gas-phase and 58 ices. Our reaction 
network contains 2479 reactions, including thermal desorption, gas-grain reactions and grain 
surface reactions. No non-thermal desorption processes (e.g. cosmic ray heating) are in- 
cluded. 

We make the simplifying assumption that radial and vertical motions of gas and dust 
are slow compared to chemical reaction timescales, and do not track the motion of gas and 
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dust through the nebula. This assumption is true for grain surface reactions and freeze- 
out/desorption, which have short timescales, but gas-phase chemistry proceeds more slowly. 
Furthermore, inward migration of cm-size particles can deplete the outer solar nebula of ice 
and concentrate each species at its condensation front (e.g. Ciesla and Cuzzi 2006). In this 
paper, we ignore relative motions of solid and gas. 

Since protostellar disks have a high optical depth to UV radiation, we neglect pho- 
toionization in the disk midplane model. However, we include cosmic rays at the standard 
interstellar flux in parts of the disk with a half-plane surface density lower than the pene- 
tration depth of 150 g cm -2 (Umebayashi and Nakano 1981). The cosmic-ray penetration 
depth is listed in Table [3] with the other free parameters in the chemical model. 

The following sections describe our treatment of gas-grain and grain surface reactions 
and the binding energies we use to calculate those reaction rates. 



4-1.1. Grain Surface Reactions 



Following Willacy et al. (2006) and Garrod and Herbst (2006), we use rate equations to 



calculate grain surface reaction rates. Rate equations work in the mean-field approximation, 
neglecting the stochastic variation of abundances on different grains. The accuracy of rate 
equations is compromised when a species reacts on the grain surface faster than it can adsorb 
or desorb: this drives the reactant abundance per grain to less than unity and the mean-field 
approximation is no longer valid (Caselli et al. 1998). However, unlike stochastic methods 
that account for the discrete nature of the grains, solving rate equations is computationally 
tractable even for large-scale simulations. For example, a model that uses direct integration 
of the master equation (the probability distribution governing reaction rates) has the number 
of equations in the simulation increasing exponentially with the number of species (Barzel 
and Biham 2007). 

The species most likely to violate the mean-field approximation is atomic hydrogen, 
which can scan grain surfaces and find reaction partners quickly due to its low atomic weight. 
While heavier atoms move across grain surfaces by thermal hopping, Hornekaer et al. (2003) 
find that H atoms often tunnel between adjacent vacancies in the grain surface lattice. This 
result conflicts with previous work by Katz et al. (1999), who found that the quantum 
tunneling mechanism is improbable. Although the question of whether quantum tunneling 
or thermal hopping is primarily responsible for hydrogen mobility on grain surfaces is not 
resolved, Ruffle and Herbst (2000) show that rate equation simulations limiting H motion 
to thermal hopping only — and thus slowing down reactions involving atomic H — show good 
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agreement with stochastic methods. We therefore assume thermal hopping regulates all grain 
surface reaction rates. 

The thermal hopping rate to an adjacent vacancy on the grain surface lattice is 

*ho P = exp (-0.3E D /kT g ) , (28) 

where Ed is the binding energy, z/ is the oscillation frequency between the ice and the grain 
surface, and T g is the grain surface temperature (assumed to be in equilibrium with the gas). 
The oscillation frequency is 

u = y/2n s E D /Tt 2 m } (29) 

where n s ~ 1.5 x 10 15 cm~ 2 is the surface density of lattice vacancies (Hasegawa et al. 1992; 
Tablepl) and m is the reactant mass. Typical oscillation frequencies are 1-3 x 10 12 s _1 (Herbst 



et al. 2005). Following Willacy et al. (2006), we assume that only atoms can move across 
grain surfaces: molecules, which have higher binding energies, are stationary. All grain 
surface reactions thus involve at least one atom. 

The grain-surface reaction rate (s _1 ) is given by 

1 _ ^hop,! T l hoPi2 , . 

rcac " Ana?n s n g ' W 

where thop.i and t hop 2 are the thermal hopping rates of reactants 1 and 2, a is the grain radius 
and n g is the number density of grains. As in the molecular cloud model, we assume a mean 
grain radius of 0.258 \xm (Table g. We take a grain number density of n g / {n-n+2nm) = 10 12 
(Table [3]), consistent with our assumption of well-mixed grains and gas. 



4-1-2. Gas-Grain Reactions 

We assume that all species can freeze out on grains except He, which has a very low 
binding energy, Ep/k = 100 K (Tielens and Hagen 1982). The freeze-out rate is 

k f = S x (Tra 2 n g ) CiV x n x , (31) 

where S x is the sticking coefficient (assumed to be 0.3 for H atoms and 1 for all other species; 
Table [3]), v x is the gas-phase velocity of species x, and n x is the gas-phase number density of 
x. The factor Cj takes into account the attraction of positively charged ions to grains which, 
on average, are negatively charged (Umebayashi and Nakano 1980): 

( Ci = 1 for neutrals 

\ C{ = 1 + e/ (akT) for single — charged positive ions. 
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In Eq. 32, k is the Boltzmann constant and e is the electron charge. Once an ion hits a 
grain, we assume it is neutralized by an adsorbed free electron. The thermal desorption rate 



[s ) at which ice mantle species are returned to the gas is 



k t = v exp (E D /kT g 



(33) 



4-1.3. Binding Energies 

The binding energies Ed for each species depend on the surface onto which the molecules 
adsorb and the mixture of ices on that surface. If each ice were a pure solid, Ed would simply 
be the sublimation energy and desorption physics could be calculated nonkinetically using 
vapor pressure. However, molecules in mixed ices may have substantially different binding 
energies than those in pure ices, so the sublimation energy approximation is inadequate. 

The most abundant astrophysical ice is H 2 0, so the binding energy of each species under 
protostellar disk conditions is primarily determined by how it behaves on/in an amorphous 



water-ice matrix. Collings et al. (2004) group binding energies into three categories based 
on their behavior when co-deposited with H2O: CO-like, water-like, and intermediate. The 
extremely volatile CO-like species include CO, N 2 , 2 , and CH 4 . Molecules of these species 
do not bond efficiently either with themselves or with water: they quickly diffuse through 
an amorphous H 2 matrix and desorb mainly at their sublimation temperatures, which are 
quite low (< 40 K), with minor amounts remaining trapped in the H 2 ice. 

Molecules of water-like species, including NH 3 and CH 3 OH, bind more strongly to water 
than to each other or to a metal substrate. These polar molecules form hydrogen bonds with 
water that keep them in the solid phase at higher temperatures than possible in their pure 



form. The data of Collings et al. (2004) show that ammonia- water mixtures are substantially 
less volatile than pure ammonia, which has a binding energy of only 1110 K (Hasegawa and 
Herbst 1993). NH 3 therefore takes on the binding energy of H 2 0, 5700 K. 

Finally, intermediate species are trapped by H 2 to some extent, but do not have the 
same mobility as CO inside a water-ice matrix. Some co-desorption with H 2 occurs, but 
most desorption takes place at the sublimation energy. These species include H 2 S, OCS and 
CH 3 CN. We have not included the partial co-desorption with water and treat each molecule 
as having only one specific desorption tmeperature. Table [4] gives the binding energy data 
used in our model. 



4.2. Preprocessing in Molecular Cloud 



The ices in protostellar disks have undergone significant chemical processing before giant 
planet formation begins. To simulate this preprocessing, we model 10 6 years of molecular 
cloud evolution to derive input abundances for the protostellar disk. The molecular cloud 
model uses the same species set and reactions as our disk model, with the addition of 
photoionization due to the interstellar UV field. 

The atoms that form the species in our simulation are H, He, C, N, O, and S. We 
also include Si adsorption, desorption and ionization as a tracer of the metal population. 
Except for hydrogen, we assume all material enters the molecular cloud in atomic form, as 
if coming from a recent supernova. We take nn/(^H + 2nH2) = 10~ 2 (Table [3]), so that 99% 
of hydrogen atoms are in H2. We assume a density of 2 x 10 4 cm -3 , a temperature of 10 K 
and an extinction Ay = 10 mag. 

Not all of the atomic CHONS material is available to form ices: some will help form the 
solid grains where ice mantles deposit. To figure out the fraction of C, N, O, and S locked 
away in grains, we use the atomic abundance ratios in Halley's comet dust grains reported 



by |Jessberger et al.| fll988[ ): C:N:0:Si:S = 814:42:890:185:72. Assuming only trace amounts 
of silicon remain in the gas phase (gas/grain < 0.001), we calculate grain fractions for C, N, 
O and S as 



\ S 1 / comet 



rWn gas =(— j I — ) (34) 



Using the solar composition of Helling et al. (2000), as in the dynamical disk model, we 



find -F grain = 0.4, 0.1, 0.2 and 0.85 for C, N, O and S respectively (Table [3J. This portion of 
each atom's inventory does not participate in chemistry and is excluded from the simulation. 
The gas and ice abundances of major species after 1 Myr of molecular cloud evolution are 
listed in Table [5] These abundances form the chemical initial conditions for our protostellar 
disk model. 

The major ice species formed in the molecular cloud are H2O, NH3, CO, and CH4. 
More than 90% of the C, N and O atoms that aren't part of refractory grains are in these 
compounds. Sulfur is sequestered in solid H 2 S after the molecular cloud phase, but this 
composition is subsequently modified during disk evolution. Minor CNO carriers are HCN, 
N2 and NO. Figure [7] shows the division of solid C, N and O atoms among the species in our 
model. 

Our molecular cloud model predicts total CO/CH 4 = 0.14. Both of these species reside 
in ice mantles. Abundances of major ices relative to H2O are listed in Table 16} Our CO/CH4 



ice ratio is consistent with the Aikawa et al. (2008) models of prestellar cores, which give 



- 22 - 



CO/CH 4 = 0.43 in an embedded core and CO/CH 4 = 0.1 in an isolated core. Likewise, 
Lodders (2003) concluded that the reactions 

3H 2 + CO -> CH 4 + H 2 

CO -> C (graphite) + C0 2 [ ' 

would deplete CO from presolar gas (note, however, that in our model CH4 forms on grains 
and not in the gas). Lodders predicts CH 4 /H 2 = 0.58, whereas we find a methane-to- water 
ratio of 0.38 (Table [6]). This discrepancy reflects the fact that we have a higher input oxygen 
abundance (to match the opacity tables calculated by Semenov [2003]). 

We note that the dominant components of interstellar ice mantles are typically H 2 0, 
NH 3 , CO, C0 2 , and CH 3 OH (Charnley and Rodgers 2008), which suggests carbon oxide 
formation in the ISM should be favored and hydrocarbon formation kinetically inhibited. 
One possible way to reduce the CH 4 /CO ratio is to lower the assumed atomic hydrogen 
abundance, which would reduce the efficiency of hydrogenated molecule formation on the 
grains at early times when atomic carbon is abundant. However, the purpose of this work is 
to determine the ice mass available for giant planet formation. Since CH 4 and CO have such 
similar condensation temperatures (41 and 35 K, respectively) and we assume O/C > 1, so 
CO cannot incorporate all available oxygen, CO/CH4 is not critical to our final result. 

We find that the most abundant nitrogen-carrying molecule is ammonia, with N 2 /NH3 = 
2.5 x 10~ 3 (both species reside in ice mantles). We calculate NH 3 /H 2 = 0.14, which matches 
the input N/O ratio. Lodders (2003) also predicted that hydrated ammonia was the major 



nitrogen reservoir in the solar nebula (see £4.3.3 for a discussion of hydrated ammonia in 



Jupiter and Saturn's feeding zones). However, Womack et al. (1992) infer N 2 /NH3 = 4 in 
the comet-forming region of the solar nebula from models of elemental nitrogen depletion in 
Comet Halley. 

Overall, our model molecular cloud favors the formation of hydrogenated C, N, and O 
compounds over oxides or diatomic molecules (except for H 2 ). We next discuss the subse- 
quent chemical evolution of the solar nebula. 



4.3. Chemical Model Results: Ice Lines 

In this section, we report the locations of important condensation fronts and the relative 
abundances of the most common ices. We use the gas and ice mixture resulting from the 
molecular cloud model as the starting point for a radial series of midplane chemical models 
that span 2 Myr of disk evolution. These models track the ice inventory of the disk midplane 
as a function of heliocentric distance and time. 
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The ices in the disk can be subdivided into two categories: (1) those that form in the 
molecular cloud and freeze out from the disk gas, and (2) those that form or are destroyed 
in the warm disk midplane. All Category 1 molecules are highly stable both as gases and 
solids: in the radiation-shielded disk midplane, their reaction probabilities are very low. The 
stability and ease of formation of these molecules makes them highly abundant, > 99% of 
the solar nebula ice mass. In Category 1 are water, methane, CO, ammonia, and HCN. In 
Category 2 are all sulfur-carrying molecules, simple hydrocarbons, N 2 , and NO. 



4-3.1. Category 1: Stable and Abundant 

Figure [8] shows time snapshots of the abundance of Category 1 ice abundances (solid 
phase only) as a function of radius. Note that the ammonia and water condensation fronts 



coincide; we will discuss the implications of this further in £ 4.3.3 The most abundant ice is 
H 2 0, followed by methane, ammonia, CO and HCN. These five ices account for 98% of the 
ice mass in regions where T < 35 K, the sublimation temperature of CO. 

In Fig. [8j we see the H2O condensation front sweeping from 6 AU at t — to 1.5 AU 
after 2 Myr. Observations of the asteroid belt place the H 2 ice line at 2.7 AU (Abe et 
al. 2000; Rivkin et al. 2002), a position it reaches at t = 5 x 10 5 yr in our model. The 
slow motion of the snow line through the inner solar system could perhaps account for the 
H 2 abundance gradient observed in the asteroid belt (Barucci et al. 1996): we find that 
icy grains are present for ~ 10 5 years longer at 3 AU than at 2 AU. The snow line crosses 
Jupiter's heliocentric distance, 5.2 AU, after only 5 x 10 4 yr, leaving ample time for the 
formation of Jupiter's core from icy planetesimals. 

If we assume Saturn formed in situ at 9.5 AU, primordial methane accretion onto both 
the giant planet and its satellites is possible. We find that grain mantles in Saturn's feeding 
zone should have CH 4 /H 2 > 0.01 after 2.3 x 10 5 yr. CH 4 is 90% solid at 9.5 AU after 
5.4 x 10 5 yr. We note that the C/H ratio in Saturn's atmosphere has recently been revised 
upward by a factor of 2 (Flasar et al. 2005). 

We find that if Uranus and Neptune formed outside 12 AU, their core-forming planetes- 
imals contained both CH 4 and CO. The composition of both planets strongly suggests that 
they accreted carbonaceous ices: (C/H)/(C/H) is 41 for Neptune and < 260 for Uranus 
(Lodders and Fegley 1994). Furthermore, Grundy et al. (2002) find that methane ice is 
widely distributed across the surface of Triton. 
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4-3.2. Category 2: Chemically Active 

In Fig. [8] we see that the typical behavior of an ice line is to move inward as the nebula 
cools. This can happen only if the molecule is relatively inert. Sulfur species are missing 
from Fig. [8] because H 2 S, the dominant sulfuric ice at t — with 99.6% of the sulfur atoms, 
does not survive in the gas phase. Likewise, while CsH n chains are stable and freezeout- 
dominated in our model, acetylene ice partially hydrogenates to form ethylene and ethane. 
Gaseous NO and N 2 participate in ion exchange reactions in the gas, reducing the abundance 
available for ice formation. In this section, we consider the sulfuric, nitric and hydrocarbon 
ice systems. 

As soon as our midplane chemical model begins, the reaction 

H 2 S + H -> HS + H 2 (36) 

starts destroying gaseous H 2 S where temperature is above ~ 50 K. The product HS then 
initiates a network of reactions that form OCS and H 2 CS. The reactions linking the main 
sulfur- carrying molecules are shown in Fig. [91 We do not find any depletion of H 2 CS once 



formed on grains. Figure 10 shows time snapshots of the abundances (solid phase only) 
of the main sulfur-carrying molecules H 2 S, OCS and H 2 CS. After 2 Myr, we see a banded 
structure where H 2 S dominates at R > 8 AU, OCS is most abundant between 5 and 8 AU 
and H 2 CS dominates at R < 5 AU. 

H 2 S remains in the outer solar nebula because the gas at R > 9 AU never gets hotter than 
the H 2 S sublimation temperature in our model. This scenario may be physically unrealistic: 
infalling molecular cloud material is likely heated in an accretion shock during disk formation. 
However, if we examine the total abundance of icy sulfur atoms, iV(H 2 S) + iV(OCS) + 
iV(H 2 CS), we see that it is conserved after OCS formation finishes at 3 x 10 5 yr. From 
Fig. [TUJ we can imagine a sulfurous ice line sweeping through the solar nebula, moving from 
8.5 AU at t = to 2.5 AU after 2 Myr. 



Figure [10] also shows the abundances of solid NO and N 2 . NO is chemically active during 
the late stages of the solar nebula, T > 10 6 yr. Even when the solar nebula temperature drops 
below the nominal sublimation temperature of NO, some thermal desorption still takes place, 



albeit at a much-reduced rate (see Eq. 33). For most species, re-freezing after desorption 
occurs faster than reactions in the gas and no ice molecules are lost. However, NO can be 
destroyed in the gas by reactions with H + . Once the half-plane surface density drops below 
the cosmic ray penetration depth at late times, the production rate of H + increases by 33% 
due to the reaction of H 2 with cosmic ray-produced photons: 

H 2 + 7 Ci^H + + H + e-. (37) 
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Thermally desorbed NO atoms are then efficiently removed from the gas before refreezing 
by the reaction sequence 

NO + H+ ->• NO+ + H 
then 

NQ+ freeze N + Q (3g) 

or 

NO + + e" ^N + O. 

Laboratory work by Oberg et al. (2005) indicates that N 2 is extremely volatile, with 
Eo/k = 790 K. This puts its sublimation temperature near the assumed ambient temper- 
ature of our protostellar disk simulation, 20 K. The smooth shape of the N 2 abundance 
curve at t = 1 yr reflects the midplane temperature's asymptotic approach to isothermality 
at 20 K. As the disk evolves, the N 2 ice line moves inward, but the overall N 2 abundance 
declines due to the gas-phase reaction 

He + + N 2 -> N+ + N + He. (39) 

The N + ion then begins a sequence of reactions with H 2 that produce NH + , N 2 H + , N 3 H + 
and finally the ammonium ion N 4 H + . Ammonium recombines with an electron to form NH 3 
+ H, so the net effect is to convert gaseous N 2 into NH3. 

Our model includes the small aliphatic hydrocarbons C 2 H 2 , C 2 H4, C 2 H6, C3H 2 and 
C3H4. We do not include longer chains in our model — we merely use these five species 
to find the location of the "oil line" of icy hydrocarbons. Lodders (2004) suggested that 
tarry compounds, rather than water, provided the bulk of the solid mass for Jupiter's core. 
This work does not truly test that hypothesis, since Lodders' scenario requires a substantial 
amount of hydrocarbons that are more refractory than water (our light oils have sublimation 
temperatures of 55-70 K). PAHs would be suitable tar candidates, but we do not include them 
in our model for lack of accurate binding energy measurements. We plan on incorporating 
PAH freezeout into future work, but we note that Jupiter falls inside the water-ice line for 
all but the first 50,000 yr of disk evolution. 



Figure 11 shows the evolution of the hydrocarbon ice lines. The most abundant hydro- 
carbons in our simulation are C3H n chains, which are at the end of our reaction sequence. 
This suggests that concatenation of small carbon chains is efficient. The small hydrocarbons 
in our model have similar sublimation temperatures, creating an "oil line" that moves from 
8 AU at t = to 3 AU after 2 Myr. 

For the C 2 H n hydrocarbons, acetylene dominates over ethylene and ethane in the trans- 
Saturnian region, R > 10 AU. The activation barrier E/k = 1210 K for the reaction 
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C2H2 + H — > C2H3 (Hasegawa et al. 1992), which has to break the strong C-C triple bond, 
means that at the cold temperatures required for acetylene to deposit on grain surfaces, 
hydrogenation proceeds slowly. Only when the disk is near the acetylene sublimation tem- 
perature (~ 55 K) does a significant amount of ethane form. Throughout most of the disk, 
ethylene is more abundant than ethane. However, an "ethylene gap" forms between 7 and 
10 AU after 3 x 10 4 yr. We explain the ethylene gap as follows: between 10 4 and 10 5 yr, 
7-10 AU region takes on mass due to an outward-moving density wave (recall the discussion 
in £3.7.1). This surface density increase causes a temperature increase of ~ 15 K, increasing 
the rate of ethylene hydrogenation. 



4-3.3. The Snow- Ammonia Line 

Since we have adapted lab results showing co-desorption of water and ammonia (Collings 
et al. 2004) for our model, we find that the ammonia and water condensation fronts coincide. 
This statement may seem controversial, and indeed, further lab experiments are needed 
to verify that hydrogen bonding between the species is strong enough to govern ammonia 
desorption. However, the composition of Jupiter, Saturn and their satellites provides strong 
evidence for primordial ammonia-rich planetesimals between 5 and 10 AU. 

Lopes et al. (2007) examined the rheological properties of a cryovolcanic flow on Titan 
with Synthetic Aperture Radar imaging and found a slurry composition most consistent 
with an ammonia-water-methanol mixture (recall that methanol is another ice found to co- 
desorb with water in lab experiments). The presence of ammonia in cryovolcanic outflows 
strongly suggests primordial ammonia accretion: late-stage volatile delivery by comets would 
provide a small surface reservoir of ammonia but not embed it in the geologically active 
subsurface. Based on the low temperatures (T < 75 K) necessary for trapping N2 in water 
ice, Owen (2000) argued that Titan's massive nitrogen atmosphere originally supplied as 



NH 3 , unmodified from its relative abundance in the outer solar nebula. 

Other Saturnian moons also show evidence of primordial ammonia accretion: Prentice 
(2007) calculates that Iapetus is 27% ammonia by mass, assuming homologous contraction 
of the Saturnian subnebula and placing Enceladus just at the stability point of liquid water. 
The source material for the rock/ice Iapetus cannot be solar nebula gas, nor could N 2 have 
provided the initial nitrogen budget: with 34% rock and 34% H2O, there would have been 
no free hydrogen reservoir for ammonia formation. Freeman et al. (2007) suggest that an 
ammonia-water ice mantle is necessary to explain the subsurface ocean driving Enceladus' 
south pole geyser. Finally, Prentice (2006) finds that the lack of compressional features on 
the surface of Rhea, which are expected to result from phase II of H 2 crystalline ice, favors 



-27- 



an ammonia-rich satellite composition (25% NH 3 by mass), since ammonia inhibits phase II 
H 2 ice formation. 

In the Jovian system, conducting subsurface oceans are required to explain the magnetic 
fields of Europa, Ganymede, and Callisto. One way to lower water's melting point and 
maintain liquid oceans is to mix in ammonia (Hussmann et al. 2006). Spohn and Schubert 
(2003) use equilibrium heat transfer models to confirm the likelihood of subsurface ammonia- 
water oceans on all Galilean satellites except Io. 

Perhaps the best piece of evidence for ammonia-containing planetesimals between 5 and 
10 AU is the nitrogen enrichment in Jupiter and Saturn's atmospheres. (N/H) / (N/H) Q = 3.3 
for Jupiter and 2-4 for Saturn (Owen and Encrenaz 2003). Such a large enrichment suggests 
a non-gaseous origin for nitrogen in both planets. Trapping N 2 or NH3 in clathrate hydrates 



is possible, but Hersant et al. (2004) find that the ammonia clathrate cannot form until 
0.9 Myr after solar nebula formation: the ammonia clathrate is more volatile than a co- 
deposited ammonia-water ice mixture. Based on the upward revision of Saturn's C/H ratio, 
Hersant et al. (2008) suggest that ammonia was present as a hydrate, rather than a clathrate, 
in Saturn's feeding zone. We recommend that the canonical "snow line" be reconceived as 
a water ice-ammonia line, and that NH3 be added to the volatile inventory of Jupiter- and 
Saturn- forming planetesimals. 



5. Solid Surface Density 

Knowing the solid-phase abundance of each ice species and the location of its conden- 
sation front, we can now calculate the solid surface density available for planet formation 
everywhere in the disk. We first calculate the gas/solid mass ratio 

G/S = (40) 

where \i is the mean molecular weight of a species, n is its abundance, and the subscripts g 
and s denote gas and solid species, respectively. We include Na, Mg, Si, Fe, Ni, Al and Ca in 
the solid inventory in solar proportions and assume these atoms form refractory compounds. 
Finally, we calculate the solid surface density as 

SsoHd = Y^ffs- (41) 



Figure 12 shows the evolution of the solid surface density distribution in the solar nebula. 
Note that as the disk evolves, the solids should grow and decouple from the gas, whereas 
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these solid surface density calculations assume grains and gas are well mixed. We see the 
large bump caused by the snow-ammonia line beginning at 6 AU and moving inward to 
1.5 AU after 2 Myr. In addition, methane adsorption creates a second local maximum in 
the solid surface density profile at 11 AU, which moves inward to 6.5 AU by the end of the 
simulation. 



Panels 1 and 2 in Fig. 12, t = 1 and 3 x 10 4 yr, show the effect of the decretion flow 
on the solid surface density profile. The solid mass between 6 and 14 AU actually increases 
during the early stages of disk evolution, at the expense of the inner 4 AU. The next panel, 
10 5 yr, shows the further solid surface density increase between 12 and 18 AU. Disk annuli 
that increase in surface density also heat up, and the methane bump consequently moves 
outward for the first 10 5 yr. 



Pollack et al. (1996) presented core accretion models of Jupiter, Saturn, and Uranus, 
beginning with solid surface densities of 10, 3, and 0.75 g cm -2 . While these simulations 
showed the plausibility of giant planet formation by core accretion, the planet formation 
timescales calculated were too long: Jupiter, Saturn and Uranus took 8, 9.6, and 16 Myr to 
form, respectively. Since runaway solid core growth rate is directly proportional to surface 
density (Eq. [T]), we see that having a high surface density speeds up giant planet formation. 

If one assumes that nearly 100% of the available solids were incorporated into planetes- 
imals, our results give starting solid surface densities that are far more conducive to giant 



planet formation than those from Pollack et al. (1996). (The £ so iid of Pollack et al. (1996) 



refers specifically to 100-km planetesimals, not smaller bodies or grains.) Table 7 lists solid 
surface density as a function of radius and time in the solar nebula. Note that these results 
assume a substantial population of small grains is always available for ice deposition. 

We consider 1.5 x 10 5 yr the extreme upper limit of the planetesimal formation timescale, 
which Hubbard and Blackman [2006] predict to be < 5 x 10 4 yr. However, we continue the 
disk simulation out to 2 Myr in order to calculate nebular gas temperature, density and 
composition, which are necessary for models of giant planet formation that build on this 
work (see, for example, the new Saturn formation simulation by Dodson- Robinson et al. 
[2008]). At Saturn's heliocentric distance, 9.5 AU, the solid surface density never drops 
below 8.7 g cm -2 during the first 1.5 x 10 5 yr of disk evolution. This large solid surface 
density should provide more than a factor of 3 speedup in Saturn's solid growth phase over 
the Pollack et al. results. 

For Uranus, the relative increase in solids is even more pronounced: at 20 AU, the 
minimum solid surface density during the first 1.5 x 10 5 yr is 3.3 g cm -2 , for a factor of 
4 speedup in runaway core growth. Moving Uranus' starting position inward to 16 AU, 
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in accordance with the Nice model of planet migration (Tsiganis et al. 2005), provides 
a starting E solid of not less than 4.9 g cm -2 , in addition to the speedup given by a higher 
angular-frequency orbit. The prospects for Neptune to form rapidly at ~ 12 AU, as predicted 
by the Nice model, are quite good: £ so iid > 7.1 g cm -2 . 



Hubickyj et al. (2005) provided an update to the Pollack et al. simulations of Jupiter's 
growth, this time assuming rapid grain settling and sublimation would lower opacities of 
the gaseous protoplanetary envelope and speed up its contraction. (Even if we assume low- 
opacity envelopes for Saturn and Uranus, they still require high solid surface densities to 
form quickly, as their solid growth stages both take > 1.5 Myr in the Pollack et al. models.) 
Hubickyj et al. found that, given a planetesimal surface density of 10 g cm -2 , Jupiter could 
reach its present mass in 2.3 Myr. Our work shows that this starting condition for Jupiter 
is not only viable but likely: we predict solid surface densities not less than 8.2 g cm -2 
with refractory material only and as high as 14.5 g cm -2 after the ice line moves through at 
5 x 10 4 yr. 



6. Conclusions 



By combining a viscously evolving protostellar disk model with a kinetic model of ice 
formation, we have calculated the time-evolving solid surface density available for giant 
planet formation. We find three results that are highly favorable for gas- and ice-giant 
formation: 



1. The total (gas+solid) surface density distribution evolves toward uniformity, with an 
outward-moving wave depositing mass in the giant planet-forming region between 5 
and 20 AU, 

2. The ammonia and water ice lines coincide, providing a mass enhancement of 7% at the 
snow line over pure water condensation, and 

3. Methane and CO condensation beyond 12 AU add mass to Uranus and potentially 
Neptune's feeding zones. 



The next step is to use a core-accretion model to calculate how fast Saturn, Uranus 
and Neptune can form with our enhanced surface densities. The same methods used by 
Hubickyj et al. (2005) should work for Saturn. However, an updated method of tracking 



planetesimal scattering may be required for feeding zones with R > 15 AU, where orbital 
velocities of planetesimals decrease to ~ 10-50 times their escape velocities and self-stirring 
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begins to affect the dynamics of the planetesimal disk. By calculating solid surface density as 
a function of radius and time self-consistently with viscous evolution of the protostellar disk, 
instead of assuming a single gas/solid ratio covers the entire giant planet-forming region of 
the solar nebula, the research presented here allows core accretion simulations to be fully 
deterministic. 

One area that requires more work is the viscous evolution of an MRI-turbulent disk 
with a dead zone, where a varies with radius. Annuli with an inactive midplane would have 
a lower column- averaged value of a, which could slow or stop the movement of mass from 
the inner solar nebula to the giant planet-forming region. However, Turner and Sano (2008) 
find that a toroidal magnetic field component produced in the midplane by shear drives a 
laminar accretion flow even through the nonturbulent dead zone. 

Since our focus is on increasing the solid mass available for giant planet cores, we have 
neglected clathrate hydrates in our chemical model, preferring to tap into the large mass 
reservoir provided by ices. However, noble gas enrichment in the giant planet atmospheres 
suggests that clathration is very important, if not for building planet cores, then at least 
for determining atmospheric composition during the late stages of accretion (Hersant et al. 
2004). The presence of PAHs in carbonaceous meteorites, cometary dust and Iapetus implies 
that these, too, might be an important solid mass reservoir. It is not yet known whether 
PAHs survive in the ISM or form primarily in the warm disk. Woods and Willacy (2007) 
present a mechanism for benzene formation in the inner 3 AU of the solar nebula. We are 
investigating the possibility of including PAH chemistry in future work. 

Finally, we note that the ALMA Design Reference Science Plan contains a proposal for 
detecting the snow line in nearby protoplanetary disks. We submit that the CO ice line 
would be an equally valuable detection, since this work strongly suggests that CO ice should 
be part of the solid inventory of the Uranus- and Neptune-building planetesimals. In a 1 
Myr-old solar nebula analog in Taurus (140 pc), the CO ice line would lie 6 resolution units 
from the star when observing the CO ( J = 3 — > 2) transition. 
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Table 1. Parameters Specifying Initial Solar Nebula 



Parameter 


Description 


Value 


Units 


M* 


star mass 


0.95 


M 


M disk 


disk mass 


0.12 


M 




inner radius 


0.3 


AU 


Rout 


outer radius 


30 


AU 




mass distribution 


E oc R- 15 


g cm -2 


a 


viscosity parameter 


0.002 


dimensionless 




ambient temperature 


20 


K 
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Table 2. Solar Composition Used for Opacity Table and Molecular Cloud Model 



Element 


Abundance* 


Solids Formed 


H 


12.00 


ice, refractory CHON 


He 


10.99 


none 


C 


8.60 


ice, refractory CHON 


N 


7.97 


ice, refractory CHON 





8.87 


rock, ice, refractory CHON 


Ne 


8.07 


none 


Mg 


7.58 


rock 


Si 


7.55 


rock 


S 


7.21 


rock, ice 


Fe 


7.51 


metal 



a Abundances are logarithmic, with log n H = 12.00 
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Table 3. Free Parameters in Chemical Model of Disk Midplane 



Parameter 


Description 


value 


units 


S cr 


cosmic ray penetration depth 


150 


g cm' 2 


n s 


surface density of reaction sites on grains 


10 15 


cm~ 2 


a 


mean grain radius 


0.258 




n a 


grain number density 


10~ 12 X n H 


cm~ 3 






J 0.3, x = H 
\ 1.0, i^H 




Sx 


sticking coefficient of species x 


dimensionless 


hh/K + 2n H 2) 


atomic hydrogen fraction 


0.01 


dimensionless 


™C,grain/™C 


fraction of carbon in grains 


0.4 


dimensionless 


"N, grain /™N 


fraction of nitrogen in grains 


0.1 


dimensionless 


™0,grain/™0 


fraction of oxygen in grains 


0.2 


dimensionless 


™S,grain/«S 


fraction of sulfur in grains 


0.85 


dimensionless 
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Table 4. Binding energies used to determine thermal desorption rates of ices 



Species 


Binding Energy (K) 


Reference 51 


H 


600 b 


CT 


C 


800 


TA 


N 


800 


TA 


O 


800 


TA 


s 


1100 


TA 


N 2 


790 


O 


NO 


1210 


HH 


CO 


1720 b 


SA 


H 2 


5700 


C 


HCN 


1760 


HH 


H 2 S 


1800 


HH 


C0 2 


2860 


SA 


OCS 


3000 


HH 


NH 3 


5700 


C 


CH 4 


1360 


HH 


H 2 CS 


2250 


HH 


C 2 H 2 


1610 


HH 



References: CT = Cazaux & Tielens (2002); 
TA = Tielens & Allamandola (1987); O = Oberg 
ct al. (2005); HH = Hasegawa & Herbst (1993); 
SA = Sandford & Allamandola (1990); C = 
Collings ct al. (2004). 

b Surface binding energy on an amorphous H 2 
ice substrate. 



Table 5. Output Abundances from Molecular Cloud Simulation 



Species Abundance 





4 


08 


'-05) a 


He 


9 


77 


'-02) 


Q 


2 


61 


'-1 11 


N 


2 


uo 


'111 


o 


3 


76 


'111 


n 


A 
I 


/I Q 


r U7) 


He + 


7 


01 


-09) 


TT+ 

H 3 





O ( 


' i nl 
r lU) 


T_T 


/i 
1 


OO 

yy 


' m *\ 
,-UlJ 


no 
L^o 


Q 
O 


u / 






2 


o2 


V -U5J 


bxOl 





UU 


v -Uo ) 


brlN u 


( 


/1 "3 

4o 




GH 2 o 


5 


27 


;-04) 


GHCN 


7 


20 


;-06) 


GHNC 


8 


74 


;-08) 


GCO2 


3 


80 


;-09) 


GNH 3 


7 


56 


(-05) 


GC 2 H 2 


1 


48 


;-07) 


GH 2 CO 


9 


80 


;-08) 


GCH4 


2 


02 


;-04) 


GC 3 H 2 


2 


94 


;-07) 


GHC3N 


1 


31 


;-08) 


GCH3OH 


1 


96 


(-09) 


GCH3CN 


4 


53 


(-09) 


GN 2 


1 


85 


(-07) 


GCH 2 CO 


9 


80 


(-10) 


GC 2 


1 


12 


(-10) 


GCH3CHO 


1 


61 


(-11) 


GC 2 H 4 


2 


38 


(-08) 


GH 2 S 


2 


42 


(-06) 


GH 2 CS 


2 


54 


(-09) 


GC 2 S 


6 


60 


(-09) 


GS0 2 


5 


06 


(-10) 


GC 3 S 


5 


19 


(-10) 


GC 2 H 6 


3 


77 


(-12) 


GC 3 H 4 


7 


14 


;-u) 



a All abundances are fractional with respect to the total number of hydrogen nuclei, rin + 2ri£j 2 . 1.00 (-10) should be read 
as 1 X 10- 10 . 

b G denotes solid-phase grain mantle species. 



-44- 



Table 6. Abundance Ratios of Major Ices After Molecular Cloud Phase 



Ice Species 


n x /n H . 2 o 


H 2 


= 1 


CH 4 


0.38 


NH 3 


0.14 


CO 


0.054 


HCN 


0.014 


NO 


0.0014 


N 2 


0.00035 
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Fig. 1. — Opacity as a function of temperature and density. The black line shows the 
unsmoothed tables of Ferguson et al. (2005) (T > 1000K)and Semenov et al. (2003) (T < 
700K). In the range 700 < T < 1000K, we take a weighted average in log-T space of the two 
tables. Finally, the green curve shows the smoothed opacity functions we use in the model. 
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Fig. 2. — Surface density as a function of radius and time in the giant planet-forming region. 
Each panel shows a snapshot of the disk surface density profile in black, with the profiles 
from previous epochs retained in gray. During the first 5 x 10 5 yr of evolution, a wave of 
local surface density enhancement propagates outward, depositing mass in the giant planet 
region of the solar nebula. 
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Fig. 3. — Surface density as a function of radius and time throughout the whole disk. Each 
panel shows a snapshot of the disk surface density profile in black, with the profiles from 
previous epochs retained in gray. The disk expands from 30 AU to > 80 AU over the course 
of the simulation. 
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Fig. 4. — Surface density as a function of time at 5, 10, 15 and 20 AU. An outward- 
propagating wave adds mass to each annulus as it passes through. The largest relative 
surface density increase, AE/ Sigmao, occurs at 10 AU, near Saturn's feeding zone. 
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Fig. 5. — Time snapshots of accretion rate as a function of radius. Negative values of M 
indicate "decretion" regions of the disk with net outward mass flow. 
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Fig. 6. — Midplane temperature as a function of radius and time in the inner 30 AU of 
the disk. Outside of 30 AU, the midplane is isothermal at 20 K, the assumed ambient 
temperature. The dashed line at 800 K shows the minimum temperature for silicate crys- 
tallization (Gail 1998); the dotted line at 160 K shows the sublimation temperature of H 2 
from the temperature-programmed desorption experiments of Collings et al. (2004); and the 
dash-dotted line at 20 K shows the assumed ambient temperature. 
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Fig. 7. — Distribution of C, N and O atoms among dust and most common ice species. 
Proportions are rounded to the nearest 1%. 
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Fig. 8. — Time snapshots of ice abundance as a function of radius for species formed mainly 
by gas freezeout. Note that the condensation fronts are not step functions: the region of 
partial ice freezeout can cover up to 3 AU. 
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Fig. 9. — Reactions linking the most abundant sulfuric ices, H 2 S, H 2 CS and OCS. The label 
"G" in front of a compound name denotes that the species is solid, adsorbed on the grains. 
The notation H2S — > HS denotes the reaction H2S + H — > HS + H2] all other reaction 
pathways can be read analogously to this example. Non-sulfuric reaction by-products, such 
as H 2 in the aforementioned reaction, are not shown. An unlabeled double arrow, such as 
the one connecting GOCS and OCS, denotes a phase change. Reaction pairs are not in 
equilibrium. 



- 55 - 



1 yr 



3 x Id 4 yr 




a 10 ia 

Radius (AU) 

10° yr 



a 10 ta 

Radius (AU) 
3 x 10 s yr 




a io ia 

Radius (AU) 
6 X lQ a yr 



a io ia 

Radius (AU) 



10' yr 




a io ia 

Radius (AU) 
1.6 x 10' yr 



a io ia 

Radius (AU) 
g x lQ* yr 




a io ia 

Radius (AU) 



a io ia 

Radius (AU) 



Fig. 10. — Abundances of solid NO, N 2 , and sulfur-carrying species as a function of radius 
and time. We find a banded sulfur ice distribution in which H 2 S dominates outside ~ 7 AU, 
OCS is most abundant between 5 and 7 AU, and H 2 CS is dominant in the inner nebula. 



- 56 - 



1 yr 



C3Hn 
■ C8H8 




r C8H6 h 

; J 






'1- 



10"* ■ 
io-* r 
10* r 

10"* 



io-' 
io-* 

10" T 
10"* 
IO"* 



3 x iff* yr 



5 10 15 

Radius (AU) 

10° yr 



E0 




10"* ■ 
io-* r 

10'* r 
10"* 




5 10 IS 

Radius (AU) 



6 x 1Q 5 yr 



20 



6 10 15 

Radius (AU) 



10° yr 



S 10 IS 

Radius (AU) 

1.6 x io° yr 



5 10 13 

Radius (AU) 

2 x lQ* yr 



5 10 IS 

Radius (AU) 



S 10 15 

Radius (AU) 



5 10 15 20 

Radius (AU) 

3 x 10 s yr 



30 





Fig. 11. — Hydrocarbon ice abundances as a function of radius and time. We group all CsH„ 
compounds together because our reaction sequence ends at C 3 H 4 , so we cannot differentiate 
between different levels of C 3 saturation. 
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Fig. 12. — Time snapshots of solid surface density available for planet formation as a function 
of radius. In each panel, the current solid surface density is plotted in black and previous 
time snapshots are retained in gray. 



